/*
  Translated from ImageJ Java Pluging into MATLAB MEX c code by
  Tom Macura <tmacura@nih.gov>
  July 21st, 2006
  Modded by Chris Coletta and Ilya Goldberg to deal with 16-bit/chanel RGB,
  and make it independent of MATLAB
*/

/* G.Landini at bham ac uk
   30/Mar/2004 released
   03/Apr/2004 resolved ROI exiting
   07/Apr/2004 added Methyl Green DAB vectors
   08/Jul/2004 shortened the code
   01/Aug/2005 added fast red/blue/DAB vectors
   02/Nov/2005 changed code to work with image stacks (DLC - dchao at fhcrc org)
   02/Nov/2005 changed field names so user-defined colours can be set within 
               macros (DLC - dchao at fhcrc org)

   This plugin implements stain separation using the colour deconvolution
   method described in:

       Ruifrok AC, Johnston DA. Quantification of histological
       staining by color deconvolution. Analytical & Quantitative
       Cytology & Histology 2001; 23: 291-299.

   The code is based on "Color separation-30", a macro for NIH Image kindly provided
   by A.C. Ruifrok. Thanks Arnout!

   The plugin assumes images generated by color subtraction (i.e. light-absorbing dyes
   such as those used in bright field histology or ink on printed paper) but the dyes
   should not be neutral grey.

   I strongly suggest to read the paper reference above to understand how to determine
   new vectors and how the whole procedure works.

   The plugin works correctly when the background is neutral (white to light grey), 
   so background subtraction and colour correction must be applied to the images before 
   processing.
 
   The plugin provides a number of "bulit in" stain vectors some of which were determined
   experimentally in our lab (marked GL), but you may have to determine your own vectors to
   provide a more accurate stain separation, depending on the stains and methods you use.
   Ideally, vector determination should be done on slides stained with only one colour
   at a time (using the "From ROI" interactive option).
 
   The plugin takes an RGB image and returns three 8-bit images. If the specimen is
   stained with a 2 colour scheme (such as H & E) the 3rd image represents the
   complimentary of the first two colours (i.e. green).
  
   Please be *very* careful about how to interpret the results of colour deconvolution
   when analysing histological images.
   Most staining methods are not stochiometric and so intensity of the chromogen
   may not correlate well with the quantity of the reactants.
   This means that intensity of the colour may not be a good indicator of
   the amount of material stained.

   Read the paper!
*/

#include <math.h>
#include <float.h> // for DBL_MAX
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <tiffio.h>
#include "readTIFF.h"



typedef struct  {
	char   *argP;
	char   label[256];
	char   isRGB;
	double MODx,MODy,MODz;
	double max,min;
} stain_t;

void usage (const char *argv0, const char *message);
void prepare_conversion_vectors(char **stain, double *q);
void getStainVecs (stain_t *stains, int nstains);
void getNormStainMODs (stain_t *stains, double *RGBwht, double *cosx, double *cosy, double *cosz);
void getStainMat(stain_t *stains, double *RGBwht, double *q);
void color_decon (float *inBuf,float *outBuf,size_t imagesize,double numbits,double *RGBwht,stain_t *stains);
void getNewMOD (float *inBuf,float *outBuf,size_t imagesize,double *RGBwht,stain_t *stains, int nstains, double *MODx1,double *MODy1,double *MODz1);


int main (int argc, char **argv) 
{
int arg,ret;
TIFF* tif;
float *inBuf,*outBuf;
void *img_stain;
uint16_t bits;
uint32_t width,height;
char outfile[256],basename[256],ext[256],error[256],*cp;
int nstains=0;
stain_t stains[3];
size_t i,j;
size_t imagesize;
uint8_t *img_stain8;
uint16_t *img_stain16;
uint32_t *img_stain32;
int minbits; // The minimum number of bits to express the maximum RGB value in the TIFF
             // Calculated by ReadRGBTIFFDataFlt
char whtLvl=0;
double MODx1[3],MODy1[3],MODz1[3];
double numbits;
double RGBwht[3];
double Rin,Gin,Bin;

	if (argc < 3) {
		usage (argv[0],"Insufficient arguments");
	}

   /*
     String [] stain={"From ROI", "H&E", "H DAB","FastRed FastBlue DAB","Methyl Green DAB", 
                      "H&E DAB", "H AEC","Azan-Mallory","Alcian blue & H","H PAS","RGB","CMY", 
                      "User values"};

    IGG: stain is specified per-channel, so -s0 H -s1 E will put the default Haem matrix 
         in MOD*[0], and the Eosin matrix in MOD*[1]. Viable stain strings are: "H", "E", 
         "DAB", "FastRed", "FastBlue", "MethylGreen", "AEC", "Azan-Mallory", "AlcianBlue", 
         "PAS", "R", "G", "B", "C", "M", "Y" To use specific user values, use 
         -s0 '123.3 456.7 234.3', -s1 '1.23 4.56 7.89'  corresponding to the RGB values.
         The parameter following s0, s1 or s2 is checked against all the stain strings, 
         and finally we try to sscanf 3 floats.  If we don't get 3 floats, then the 
         parameter is malformed.
   */
  
	// Get the CLI params
	for (arg = 1; arg < argc-1; arg++) {
		// Stain specification (will be parsed/checked later)
		if ( !strcmp (argv[arg],"-s")  || !strcmp (argv[arg],"-sOD")) {
			nstains++;
			if (nstains < 4) {
				stains[nstains-1].argP = argv[arg+1];
			// Determine if an RGB triplet is an OD representation of stains
				if ( sscanf (argv[arg+1],"%lf,%lf,%lf",&Rin,&Gin,&Bin) == 3  && !strcmp (argv[arg],"-s"))
					stains[nstains-1].isRGB = 1;
				else 
					stains[nstains-1].isRGB = 0;
			} else {
				usage (argv[0],"A maximum of 3 stains are supported");
			}
		// Get the white-level
		} else if ( !strcmp (argv[arg],"-w") ) {
			if ( sscanf (argv[arg+1],"%lf,%lf,%lf",&Rin,&Gin,&Bin) != 3 ) {
				usage (argv[0],"-w parameter is malformed. Three white levels expected (R,G,B)");
			} else {
				whtLvl = 1;
				RGBwht[0] = Rin;
				RGBwht[1] = Gin;
				RGBwht[2] = Bin;
				fprintf (stdout,"white-level = (%lf, %lf, %lf)\n",RGBwht[0],RGBwht[1],RGBwht[2]);
			}
		}
	}
 
	if (nstains < 2) usage (argv[0],"A minimum of 2 stains are supported");
	getStainVecs (stains, nstains);

	
  // suppress warnings
  TIFFSetWarningHandler (NULL);
  for (arg = 1; arg < argc; arg++) 
  {
   	// skip any command-line options and their parameters.
  	if (*(argv[arg]) == '-') {
  		arg++;
  		continue;
  	}

    fprintf (stdout,"File: %s\n",argv[arg]);
    tif = TIFFOpen(argv[arg],"r");
    if (!tif) {
      fprintf (stderr,"Could not open %s\n",argv[arg]);
      continue;
    }
    
    // Check if its proper RGB
    if ( (ret = isRGB (tif)) < 0) {
      GetReadTIFFError (ret,error);
      fprintf (stderr, "Malformed RGB TIFF %s: %s\n",argv[arg],error);
      TIFFClose (tif);
      continue;
    } else if (ret == 0) {
      fprintf (stderr, "TIFF %s is Grayscale\n",argv[arg]);
      TIFFClose (tif);
      continue;
    }
    
    TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &bits);
    TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &width);
    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &height);

    imagesize = width * height;

    /* Create an array of floats in which to put the pixel's RGB values.
       A 16-bit image will contain values from 0 to 65 535 */
    if ( (inBuf = (float *) malloc( imagesize * 3 *sizeof(float) )) == NULL) {
		fprintf (stderr,"Error allocating memory for file %s (inBuf)\n",argv[arg]);
		TIFFClose (tif);
		continue;
    }
    // allocate outbut buffer also floats (one for each of 3 stains)
    if ( (outBuf = (float *) malloc( imagesize * 3 *sizeof(float) )) == NULL) {
		fprintf (stderr,"Error allocating memory for file %s (outBuf)\n",argv[arg]);
		free (inBuf);
		TIFFClose (tif);
		continue;
    }
    // out image is one stain at a time, so we only need one.
	if ( (img_stain = malloc( imagesize * (bits/8) )) == NULL) {
		fprintf (stderr,"Error allocating memory for file %s (imgStain)\n",argv[arg]);
		free (inBuf);
		free (outBuf);
		TIFFClose (tif);
		continue;
	}

	// Read the tiff into the float buffer
    if ( (ret = ReadRGBTIFFDataFlt(tif,inBuf,&minbits)) != 0) {
		fprintf (stderr,"Error reading file %s. Return: %d\n",argv[arg],ret);
		free (inBuf);
		free (outBuf);
		free (img_stain);
		TIFFClose (tif);
		continue;
    }

	TIFFClose (tif);

	numbits = pow(2, bits) - 1.0;
	// If white level was not specified, set it to numbits.
	if (!whtLvl) RGBwht[0] = RGBwht[1] = RGBwht[2] = numbits;

	// deconvolve the colors into outBuf
	color_decon (inBuf, outBuf, imagesize, numbits, RGBwht, stains);

	// Get new MOD vectors
	getNewMOD (inBuf, outBuf, imagesize, RGBwht, stains, nstains, MODx1, MODy1, MODz1);


	// If only two stains specified, a third image will be created the remainder
	// when the two stains are deconvolved out of the source image.
	// Create an informative "stain name" for the third image to be used in its
	// file name
	if( nstains == 2 )
		sprintf( stains[2].label, "remainder_%s%s", stains[0].label, stains[1].label );
	


    
    // Get the tiff's basename and extension
    strncpy (basename,argv[arg],255);
    if ( (cp = strrchr(basename,'.')) != NULL) {
		*cp++ = '\0';
		strncpy (ext,cp,255);
	}

    // scale and write the deconvolved images.
		for (i = 0; i < 3; i++) {
			float *pixP = outBuf+(i*imagesize);
			double output;
			// Always three files output.
			snprintf (outfile,255,"%s_%s.%s",basename,stains[i].label,ext);

			// Three possible pixel types.
			img_stain8  = img_stain;
			img_stain16 = img_stain;
			img_stain32 = img_stain;

			// Scale the result depending on bit depth.
			for (j=0; j<imagesize;j++){
			// Scale the result based on max stain value
// 				output = ( (*pixP++) / stains[i].max ) * numbits;
			// Scale the result to the max possible pixel value
				output = (*pixP++) * numbits;
			// The range of the stain is scaled to numbits for defined stains
// 				if (stains[i].argP) output = ( ((*pixP++) - stains[i].min) / (stains[i].max - stains[i].min) ) * numbits;
			// Scale the result to the max possible pixel value for undefined stains
//				else output = (*pixP++) * numbits;
			// At this point, scaled or not, the result is clipped.
				if (output > numbits) output = numbits;
				if (output < 0.0) output = 0.0;
				switch (bits) {
					case 8:
					*img_stain8++ = floor(output+.5);
					break;
	
					case 16:
					*img_stain16++ = floor(output+.5);
					break;
	
					case 32:
					*img_stain32++ = floor(output+.5);
					break;
				}
			}

			// Write out the stain
			if ( (ret = WriteTIFFData (outfile,img_stain,bits,height, width)) != 0) {
				fprintf (stderr,"Error writing file %s. Return: %d\n",outfile,ret);
				free (inBuf);
				free (outBuf);
				free (img_stain);
				continue;
			}
		}
		free (inBuf);
		free (outBuf);
		free (img_stain);

	} // end foreach arc argv
	return 0;
}


void color_decon (float *inBuf,float *outBuf,size_t imagesize,double numbits,double *RGBwht,stain_t *stains) {
	float *Rp, *Gp, *Bp;
	double R, G, B, Rlog, Glog, Blog, Rscaled, Gscaled, Bscaled, output;
	size_t i,j;
	double max[3],min[3];
	double q[9];

	Rp = inBuf;
	Gp = inBuf+imagesize;
	Bp = inBuf+(imagesize*2);
    
    max[0] = max[1] = max[2] = 0.0;
    min[0] = min[1] = min[2] = DBL_MAX;
    
    // Get the deconvolution matrix
	getStainMat(stains, RGBwht, q);

	for (i=0; i<imagesize;i++){
	/* log transform the RGB data */
		R = (double) *Rp++;
		G = (double) *Gp++;
		B = (double) *Bp++;
// From ImageJ Plug-in.  We're not scaling yet.
// 		Rlog = -((numbits*log((R+1)/RGBwht[0]))/logNumbits);
// 		Glog = -((numbits*log((G+1)/RGBwht[1]))/logNumbits);
// 		Blog = -((numbits*log((B+1)/RGBwht[2]))/logNumbits);

// Avoid taking a log of 0 when input value is 0.
// Originally, 1 wasn't added to the white-level making maximum input value come out non-zero
// 		Rlog = -log((R+1)/(RGBwht[0]));
// 		Glog = -log((G+1)/(RGBwht[1]));
// 		Blog = -log((B+1)/(RGBwht[2]));
// Adding 1 to both I and I0 makes the logs better
		Rlog = -log((R+1)/(RGBwht[0]+1));
		Glog = -log((G+1)/(RGBwht[1]+1));
		Blog = -log((B+1)/(RGBwht[2]+1));

		for (j=0; j<3; j++){
		/* Get the relative stain concentrations */
			Rscaled = Rlog * q[j*3];
			Gscaled = Glog * q[j*3+1];
			Bscaled = Blog * q[j*3+2];
// From ImageJ plug-in.  We're not scaling the output yet.
//			output = exp(-((Rscaled + Gscaled + Bscaled) - 255.0) * logNumbits / numbits);
			output = exp( -(Rscaled + Gscaled + Bscaled) );
// From ImageJ plug-in.  We're not rounding the output yet.
//			outInt=(int)floor(output+.5);

        	*(outBuf+i+(j*imagesize)) = (float)output;
			if (output > max[j]) max[j] = output;
			if (output < min[j]) min[j] = output;
//DEBUG
// if (i < 10) {
// fprintf (stderr,"%9.4lf",output);
// if (j == 2) fprintf (stderr,"\n");
// else fprintf (stderr,",");
// }
      }
    } // end for( imagesize )
    
    // Set the stains max+min
	stains[0].max = max[0];
	stains[1].max = max[1];
	stains[2].max = max[2];
	stains[0].min = min[0];
	stains[1].min = min[1];
	stains[2].min = min[2];
}


void getNewMOD (float *inBuf,float *outBuf,size_t imagesize,double *RGBwht, stain_t *stains,int nstains,double *MODx1,double *MODy1,double *MODz1) {
size_t i,j;
double min, max;
double MODx0[3],MODy0[3],MODz0[3];

	getNormStainMODs (stains, RGBwht, MODx0, MODy0, MODz0);

	// Iterate to minimize distance between MOD0 and MOD1
	// Calculate new MOD1
	for (i = 0; i < nstains; i++) {
		// Skip undefined stains.
		if (!stains[i].argP) continue;
		// Set threshold for considering new RGBs
		// N.B.: high-stain = low pixel value
		float thresh,*pixP,pix,R,G,B;
		unsigned long nPix;

		nPix=0;
		thresh = stains[i].min + (0.1 * (float)(stains[i].max - stains[i].min));
		pixP = outBuf+(i*imagesize);
		R = G = B = 0.0;
		min = stains[i].min;
		max = stains[i].max;
//fprintf (stderr,"%s min=%lf, max=%lf thresh: %9.4f\n", stains[i].label, min, max, thresh);

		for (j=0; j<imagesize;j++){
			pix = *pixP++;
//if (j < 10) {
//fprintf (stderr,"pix %lu: %9.4f\n",(unsigned long)j,pix);
//}
			if (pix < thresh) {
				float weight = 1.0 - ( (pix - min) / (max - min) );
				nPix++;
				R += (*(inBuf+j)) * weight;
				G += (*(inBuf+j+imagesize)) * weight;
				B += (*(inBuf+j+(imagesize*2))) * weight;
//if (nPix < 10) {
//fprintf (stderr,"%s pix: %9.4f weight: %9.4f\n",stains[i].label,pix,weight);
//}
			}
		}
		R /= (float)nPix;
		G /= (float)nPix;
		B /= (float)nPix;
		MODx1[i] = -log((R+1)/(RGBwht[0]+1));
		MODy1[i] = -log((G+1)/(RGBwht[1]+1));
		MODz1[i] = -log((B+1)/(RGBwht[2]+1));
		double leng = sqrt ( (MODx1[i]*MODx1[i]) + (MODy1[i]*MODy1[i]) + (MODz1[i]*MODz1[i]) );
		MODx1[i] /= leng;
		MODy1[i] /= leng;
		MODz1[i] /= leng;

		float Rd = MODx1[i] - MODx0[i];
		float Gd = MODy1[i] - MODy0[i];
		float Bd = MODz1[i] - MODz0[i];
		fprintf (stdout,"%s New RGBs: %lf,%lf,%lf (ODs: %lf,%lf,%lf -> %lf,%lf,%lf) R=%lf\n",stains[i].label,
			R, G, B,
			MODx0[i],MODy0[i],MODz0[i],
			MODx1[i],MODy1[i],MODz1[i],
			sqrt ( (Rd*Rd) + (Gd*Gd) + (Bd*Bd) ) );
	}
}

void getNormStainMODs (stain_t *stains, double *RGBwht, double *cosx, double *cosy, double *cosz) {
int i;
double leng;

	for (i=0; i<3; i++){
  		
	  	// set the cos based on MOD, accounting for wether MOD is RGB values or ODs
  		if (stains[i].isRGB) {
			cosx[i] = -log((stains[i].MODx+1)/(RGBwht[0]+1));
			cosy[i] = -log((stains[i].MODy+1)/(RGBwht[1]+1));
			cosz[i] = -log((stains[i].MODz+1)/(RGBwht[2]+1));
  		} else {
			cosx[i] = stains[i].MODx;
			cosy[i] = stains[i].MODy;
			cosz[i] = stains[i].MODz;
  		}
	  	// Normalize vector length
		leng=sqrt(cosx[i]*cosx[i] + cosy[i]*cosy[i] + cosz[i]*cosz[i]);
		if (leng >  DBL_MIN){
			cosx[i] /= leng;
			cosy[i] /= leng;
			cosz[i] /= leng;
		} else {
			cosx[i]= cosy[i] = cosz[i] = 0.0;
		}
	}
}

void getStainMat(stain_t *stains, double *RGBwht, double *q) {
double leng, A, V, C;
double cosx[3];
double cosy[3];
double cosz[3];

	getNormStainMODs (stains, RGBwht, cosx, cosy, cosz);

  /* translation matrix */
  if (cosx[1]==0.0){ /* 2nd colour is unspecified */
    if (cosy[1]==0.0){
      if (cosz[1]==0.0){
        cosx[1]=cosz[0];
        cosy[1]=cosx[0];
        cosz[1]=cosy[0];
      }  
    }
  }

  if (cosx[2]==0.0){ /* 3rd colour is unspecified - calculate complementary color*/
    if (cosy[2]==0.0){
      if (cosz[2]==0.0){
        if ((cosx[0]*cosx[0] + cosx[1]*cosx[1])> 1)
          cosx[2]=0.0;
        else
          cosx[2]=sqrt(1.0-(cosx[0]*cosx[0])-(cosx[1]*cosx[1]));

        if ((cosy[0]*cosy[0] + cosy[1]*cosy[1])> 1)
          cosy[2]=0.0;
        else
          cosy[2]=sqrt(1.0-(cosy[0]*cosy[0])-(cosy[1]*cosy[1]));

        if ((cosz[0]*cosz[0] + cosz[1]*cosz[1])> 1)
          cosz[2]=0.0;
        else
          cosz[2]=sqrt(1.0-(cosz[0]*cosz[0])-(cosz[1]*cosz[1]));
      }
    }
  }

  leng= sqrt(cosx[2]*cosx[2] + cosy[2]*cosy[2] + cosz[2]*cosz[2]);

  cosx[2]= cosx[2]/leng;
  cosy[2]= cosy[2]/leng;
  cosz[2]= cosz[2]/leng;

  /* matrix inversion */
  A = cosy[1] - cosx[1] * cosy[0] / cosx[0];
  V = cosz[1] - cosx[1] * cosz[0] / cosx[0];
  C = cosz[2] - cosy[2] * V/A + cosx[2] * (V/A * cosy[0] / cosx[0] - cosz[0] / cosx[0]);
  q[2] = (-cosx[2] / cosx[0] - cosx[2] / A * cosx[1] / cosx[0] * cosy[0] / cosx[0] + cosy[2] / A * cosx[1] / cosx[0]) / C;
  q[1] = -q[2] * V / A - cosx[1] / (cosx[0] * A);
  q[0] = 1.0 / cosx[0] - q[1] * cosy[0] / cosx[0] - q[2] * cosz[0] / cosx[0];
  q[5] = (-cosy[2] / A + cosx[2] / A * cosy[0] / cosx[0]) / C;
  q[4] = -q[5] * V / A + 1.0 / A;
  q[3] = -q[4] * cosy[0] / cosx[0] - q[5] * cosz[0] / cosx[0];
  q[8] = 1.0 / C;
  q[7] = -q[8] * V / A;
  q[6] = -q[7] * cosy[0] / cosx[0] - q[8] * cosz[0] / cosx[0];
}


void getStainVecs (stain_t *stains, int nstains) {
int i;
char errmsg[256],*labP,stainArg[256];
double Rin,Gin,Bin;

	for (i = 0; i < nstains; i++) {
		if (stains[i].argP)
			snprintf (stainArg,256,"%s",stains[i].argP);
		else
			*stainArg = '\0';

		// pull out the label if there is one
		if ( (labP = strrchr (stainArg,':')) ) {
			*labP++ = '\0';
			snprintf (stains[i].label,256,"%s",labP);
		} else {
			*(stains[i].label) = '\0';
		}

		if (!(*stainArg)) {
			/*  Zero matrix */
			stains[i].MODx = 0.0;
			stains[i].MODy = 0.0;
			stains[i].MODz = 0.0;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"%d",i);
		} else if (!strcmp(stainArg,"H")) {
			/* GL Haem matrix */
			stains[i].MODx = 0.644211; /* 0.650; */
			stains[i].MODy = 0.716556; /* 0.704; */
			stains[i].MODz = 0.266844; /* 0.286; */
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"H");
  		} else if (!strcmp(stainArg,"E")) {
			/* GL Eos matrix */
			stains[i].MODx = 0.092789; /* 0.072; */
			stains[i].MODy = 0.954111; /* 0.990; */
			stains[i].MODz = 0.283111; /* 0.105; */
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"E");
  		} else if (!strcmp(stainArg,"DAB")) {
			/* DAB matrix */
			stains[i].MODx = 0.268;
			stains[i].MODy = 0.570;
			stains[i].MODz = 0.776;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"D");
  		} else if (!strcmp(stainArg,"FastRed")) {
			/* fast red */
			stains[i].MODx = 0.21393921;
			stains[i].MODy = 0.85112669;
			stains[i].MODz = 0.47794022;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"FR");
  		} else if (!strcmp(stainArg,"FastBlue")) {
			/* fast blue */
			stains[i].MODx = 0.74890292;
			stains[i].MODy = 0.60624161;
			stains[i].MODz = 0.26731082;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"FB");
  		} else if (!strcmp(stainArg,"MethylGreen")) {
			/* MG matrix (GL) */
			stains[i].MODx = 0.98003;
			stains[i].MODy = 0.144316;
			stains[i].MODz = 0.133146;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"MG");
  		} else if (!strcmp(stainArg,"AEC")) {
			/* AEC matrix */
			stains[i].MODx = 0.2743;
			stains[i].MODy = 0.6796;
			stains[i].MODz = 0.6803;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"AEC");
  		} else if (!strcmp(stainArg,"AnilineBlue")) {
			/* Aniline Blue (AZAN) */
			stains[i].MODx = .853033;
			stains[i].MODy = .508733;
			stains[i].MODz = .112656;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"AnB");
  		} else if (!strcmp(stainArg,"Azocarmine")) {
 			/* Azocarmine matrix */
			stains[i].MODx = 0.070933;
			stains[i].MODy = 0.977311;
			stains[i].MODz = 0.198067;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"Ac");
  		} else if (!strcmp(stainArg,"AlcianBlue")) {
			/* GL Alcian Blue matrix */
			stains[i].MODx = 0.874622;
			stains[i].MODy = 0.457711;
			stains[i].MODz = 0.158256;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"AlB");
  		} else if (!strcmp(stainArg,"PAS")) {
			/* GL PAS matrix */
			stains[i].MODx = 0.175411;
			stains[i].MODy = 0.972178;
			stains[i].MODz = 0.154589;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"P");
  		} else if (!strcmp(stainArg,"R")) {
			/* R */
			stains[i].MODx = 0.0;
			stains[i].MODy = 1.0;
			stains[i].MODz = 1.0;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"R");
  		} else if (!strcmp(stainArg,"G")) {
			/* G */
			stains[i].MODx = 1.0;
			stains[i].MODy = 0.0;
			stains[i].MODz = 1.0;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"G");
  		} else if (!strcmp(stainArg,"B")) {
			/* B */
			stains[i].MODx = 1.0;
			stains[i].MODy = 1.0;
			stains[i].MODz = 0.0;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"B");
  		} else if (!strcmp(stainArg,"C")) {
			/* C */
			stains[i].MODx = 1.0;
			stains[i].MODy = 0.0;
			stains[i].MODz = 0.0;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"C");
  		} else if (!strcmp(stainArg,"M")) {
			/* M */
			stains[i].MODx = 0.0;
			stains[i].MODy = 1.0;
			stains[i].MODz = 0.0;
 			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"M");
 		} else if (!strcmp(stainArg,"Y")) {
			/* Y */
			stains[i].MODx = 0.0;
			stains[i].MODy = 0.0;
			stains[i].MODz = 1.0;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"Y");
  		} else if ( sscanf (stainArg,"%lf,%lf,%lf",&Rin,&Gin,&Bin) == 3 ) {
  			stains[i].MODx = Rin;
  			stains[i].MODy = Gin;
  			stains[i].MODz = Bin;
			if (! *(stains[i].label) ) snprintf (stains[i].label,256,"%d",i);
  		} else {
			sprintf( errmsg,"Unrecognized stain %s \n",stainArg); 
  			usage (NULL, (const char*)errmsg );
  		}  		
	}
}

void usage (const char *argv0, const char *message) {
	if (message && *message) fprintf (stderr,"%s\n",message);
	fprintf (stderr,"Usage:\n");
	if (argv0) fprintf (stderr,"%s",argv0);
	else fprintf (stderr,"colorDecon");
	
	fprintf (stderr," -s <stain>[:<label>] [-s <stain>[:<label>]] [-s <stain>[:<label>]] [-w n,n,n]\n"
					"  <stain> is one of 'H', 'E', 'DAB', 'FastRed', 'FastBlue', 'MethylGreen', 'AnilineBlue','Azocarmine',\n"
					"    'AlcianBlue', 'PAS', 'R', 'G', 'B', 'C', 'M', 'Y'\n"
					"    OR, a set of RGB values separated by commas, for example:\n"
					"      -s 1.23,4.56,7.89 -s 123,456,789:St1 -s H:Hem\n"
					"      -sOD may be used instead of -s to specify stains in optical densities instead of RGB pixel values\n"
					"    If only 2 stains are specified, the third file (file_2) contains the 'remainder'\n"
					"    A maximum of three stains can be specified\n"
					"  <label> will be appended to the stain's output file base-name after a '_'.\n"
					"    Default labels will be used if not specified.\n"
					"    The above example with a 'file.tif' filename will result in:\n"
					"    file_0.tiff, file_St1.tif and file_Hem.tif\n"
					"  -w Set white-level by specifying RGB values separated by commas\n"
	);
	exit (-1);	
}
